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OCCURRENCE  AND  POTENTIAL  IMPORTANCE  OF  SALINE  PERMAFROST  IN  ALASKA 


T.  E.  Osterkamp 
Geophysical  Institute 
University  of  Alaska 
Fairbanks,  Alaska 


Introduction 

Saline  permafrost  is  defined  to  be  permafrost  with  part  or  all  of  its  water 
content  unfrozen  because  of  the  salinity  of  the  pore  water.  However,  all  soil  pore 
water  contains  ions  thus  making  all  permafrost  saline.  Consequently,  there  is  a  need 
for  a  working  definition  of  saline  permafrost  for  geotechnical  purposes.  In  this 
paper,  permafrost  with  a  pore  water  salinity  greater  than  1  ppt  (part  per  thousand) 
or  an  electrical  conductivity  greater  than,  170  mS  m  ^  will  be  considered  saline. 

This  paper  addresses  the  occurrence  of  saline  permafrost,  primarily  in  Alaska, 
its  role  in  a  few  scientific  and  engineering  problems,  and  suggests  that  saline 
permafrost  may  be  much  more  common  and  more  impoitant  than  previously 
thought. 

Subsea  Permafrost 

Subsea  permafrost  is  the  most  obvious  and  perhaps  most  extensive  formation 
of  saline  permafrost.  It  has  been  found  in  the  submerged  continental  shelves  of  the 
Arctic  and  Antarctic  land  masses  where  pore  water  salinities  of  shelf  sediments  may 
exceed  that  of  the  overlying  seawater,  especially  in  shallow  water  .  These  increased 
salinities  are  a  result  of  the  infiltration  of  salt  derived  from  (rejected  by)  growing  sea 
ice  (Osterkamp  et  al.,  1989).  The  study  of  subsea  permafrost  dynamics  (i.e.,  its 
response  to  repetitive  periods  of  submergence  and  emergence)  is  a  scientific 
problem  in  its  own  right.  Development  of  offshore  petroleum  reserves  will  require 
an  increased  knowledge  of  the  geotechnical  properties  of  subsea  permafrost  and  of 
heat  and  matter  flow  processes  in  it. 

Coastal  Regions 

Saline  permafrost  might  be  expected  in  coastal  areas  as  a  result  of  airborne 
salts,  seawater  flooding,  and  variations  in  ground  surface  and  sea  levels  over 
geologic  time.  Such  occurrences  have  been  found  at  Nome,  Kotzebue,  Barrow, 
Prudhoe  Bay  and  at  other  sites  along  Alaska's  extensive  coast.  Extremely  high  pore 
water  salinities  have  been  found  at  the  surface,  within,  and  under  both  continuous 


and  discontinuous  permafrost  in  coastal  areas  (Osterkamp  and  Payne.  1981). 
Salinities  exceeding  20  ppt  have  been  found  within  20  m  of  the  surface  (Figure  1). 
For  comparison,  annual  sea  ice  typically  has  salinities  of  5  to  10  ppt  which  gives  some 
idea  of  the  potential  effects  on  the  geotechnical  and  physical  properties  of  the 
permafrost. 

Interior  Regions 

Observations  of  saline  permafrost  in  the  Alaska  Interior  have  been 
accumulating.  Holes  drilled  at  Fairbanks.  Eagie.  Glennallen.  and  other  sites  along 
the  pipeline  corridor  have  found  that  the  electrical  conductivity  of  the  pore  water  iii 
the  permafrost  typically  ranges  from  1 00  to  300  mS  m-i .  Some  sites  have  been  found 
with  values  in  the  500  to  1000  mS  m-i  range  (Figure  2).  At  Happy  Valley  in  the 
foothills  north  of  the  Brooks  Range,  saline  permafrost  (conductivity  2.2  x  1 03  mS  m-i) 
was  found  under  an  =  1 2m  thick  layer  of  massive  ice.  There  have  also  been  reports 
of  saline  permafrost  in  the  cold,  arid  regions  of  China.  Soviet  Union,  and  Antarctica. 
The  papers  presented  at  this  workshop  show  that  it  is  relatively  common  in  Canada. 

Potential  Importance 

Why  be  concerned  about  the  effects  of  salty  pore  water  in  permafrost?  Salts 
modify  the  chemical  potential  of  the  pore  water  in  the  permafrost.  As  a  result,  there 
are  changes  in  the  phase  equilibria  conditions,  especially  freezing  point,  water  film 
characteristics,  and  volume  of  unfrozen  pore  water  at  sub-zero  temperatures.  The 
salts  may  also  modify  the  soil  and  ice  structure;  that  of  both  pore  ice  and  ice  lenses. 
Salts  can  also  be  expected  to  be  mobile,  moving  or  concentrating  under  the 
influence  of  moisture,  chemical,  temperature,  and  hydraulic  gradients.  While  the 
details  of  these  effects  are  poorly  known,  it  is  known  that  saline  pore  water  in 
permafrost  causes  marked  changes  in  its  physical  and  mechanical  properties  and  in 
heat  and  matter  flow  processes  especially  during  freezing  or  thawing. 

Considering  the  widespread  occurrence  of  saline  permafrost,  we  can  expect  to 
encounter  a  growing  number  of  scientific,  engineering,  environmental,  and 
agricultural  problems  that  require  an  understanding  of  saline  permafrost  for  their 
rational  assessment  and  evaluation. 

Summary 

Evidence  from  subsea  permafrost,  permafrost  in  coastal  regions  of  the  Arctic 
and  Antarctic,  and  in  the  cold  and  sometimes  arid  regions  of  China,  the  Soviet  Union. 


Antarctica,  and  North  America,  indicates  that  saline  permafrost  is  very  widespread  - 
perhaps  more  common  than  non-saline  permafrost!  This  salty  permafrost  is  found 
near  the  surface  at  depths  of  geotechnical  interest,  within,  and  under  continuous 
and  discontinuous  permafrost.  We  know  that  these  salts  can  produce  significant 
changes  in  permafrost  properties  and  processes  and  that  it  is  possible  for  these  salts 
to  become  mobile  under  potential  energy  gradients.  It  is  concluded  that  saline 
permafrost  will  become  increasingly  important  in  the  rational  assessment  of 
scientific,  engineering,  environmental,  and  agricultural  problems  of  the  future. 
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Figure  1.  Electrical  conductivity  profile  for  pore  water  in 
permafrost  near  the  West  Dock,  Prudhoe  Bay,  Alaska. 
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Figure  2.  Electrical  conductivity  profile  for  pore  water  in 
permafrost  near  Livengood,  Alaska. 
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A  Model  for  the  Thermal  Regime  of  Permafrost 
Within  the  Depth  of  Annual  Temperature  Variations 


T.  Zhang,  T.  Ostarkamp,  and  J.  P.  Gosink 
Geophysical  Institute 
University  of  Alaska  Fairbanks 
Fairbanks,  Alaska 


ABSTRACT 

Temperature  changes  in  permafrost  within  the  depth  of  annual 
temperature  variations  are  often  modeled  as  a  pure  heat 
conduction  problem  in  a  homogeneous  semi-infinite  region  with  an 
equilibrium  initial  condition  and  a  harmonic  furction  of  time 
representing  seasonal  changes  in  subsurface  temperature. 
Continuous  temperature  measurements  at  the  permafrost  table  show 
that  the  temperature  there  can  be  better  approximated  by  a  step 
or  truncated  cosine  function.  Analytical  solutions  of  this  heat 
conduction  problem  have  been  obtained  for  both  step  and  truncated 
cosine  boundary  conditions.  A  temperature  oscillation  in  the 
form  of  a  step  or  truncated  cosine  function  at  the  permafrost 
table  gradually  becomes  sinusoidal  with  depth. 

INTRODUCTION 

The  thermal  regime  of  permafrost  within  the  depth  of  annual 
temperature  variations  is  of  importance  for  engineering  design  in 
cold  regions.  Soil  temperatures  are  often  modeled  as  a  heat 
conduction  problem  in  a  semi-infinite  region  with  a  periodic 
surface  temperature  due  to  the  periodic  heating  (both  daily  and 
annually) .  Temperature  variations  at  the  ground  surface  are 
usually  approximated  as  a  sinusoidal  function  of  time  and 
homogeneous  soil  materials  are  assumed  (Carslaw  and  Jaeger,  1959; 
Lunardini,  1981)^^'^^.  The  solutions  are  applied  to  investigate 
daily  and  annual  soil  ten^erature  variations  {Oke,  1987)^^^,  and 
permafrost  problems  in  cold  regions  (Lunardini,  1981)  . 
However,  for  a  wet  active  layer,  the  thermal  history  at  the 
permafrost  table  cannot  be  exactly  described  by  a  sinusoidal 
boundary  condition  since  the  temperature  there  is  constrained  at 
or  near  0*C  for  most  of  the  thaw  period  (Lachenbruch  et  al., 
1962;  Zhang,  1989)^^' 
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Temperatures  in  permafrost  are  directly  related  to  the 
temperature  variations  at  the  permafrost  table  and  this  supports  the 
choice  of  the  permafrost  table  as  the  upper  boundary  to  investigate 
the  thermal  regime  of  permafrost.  Measurements  show  that  the 
temperature  variations  at  the  permafrost  table  can  be  approximated 
by  a  step  or  truncated  cosine  function  (Lachenbruch  et  al.,  1962; 
Zhang,  1989)^^'^^.  This  paper  summarizes  analytical  solutions  to 
the  one-dimensional  heat  conduction  equation  in  permafrost  when  the 
temperature  variations  at  the  permafrost  table  are  represented  by  a 
step  or  a  truncated  cosine  function. 


MATHEMATICAL  MODELS 


A  step  function,  representation  of  the  temperature  at  the 

permafrost  table  is,  for  the  I'th  cycle. 


^(0 


0  lP<t<IP+Fi 
To  lP+Pi<t<lP  +  P2 
0  lP+P2<t<iM)P 


(1) 


where  t  is  the  time  in  days,  P  is  the  period  (one  year) ,  and  and 
are  the  turn  points  in  the  step  function.  Equation  (1)  can  be 
expanded  in  a  Fourier's  series 


«• 

m  -  +  X  A^cos  {if  (2/ .  /»2  •  Pi)] 

^  •ml  '  ^  (2) 

where  Q)*2x/P  and 

=  ,3, 

The  first  term  on  the  right  hand  side  of  (2)  is  the  mean  annual 
temperature  at  the  permafrost  table,  T.  Then,  the  mathematical 
model  can  be  written  in  the  form 


df  3x2 


(4) 


T^x,  O)*!^  G„x 


(5) 


•ml  (6) 
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(7) 


?L^Go 

X-*oodx 


where  D  is  the  thermal  diffusivity  in  units  of  m^yr~^,  x  is  the 
depth  in  m  and  Go  is  the  geothermal  gradient  in  ®C 

The  solution  for  (4)  to  (7)  can  be  written  in  the  form 


T  (i,t)  =  7-  +  Gol  +  X  ('^2t  -  Pi  ■  P,)  -  dC 


(8) 

wnere  K  «  inci/2D.  The  last  term  on  the  right  hand  side  of  (8)  is  the 
transient  disturbance  caused  by  starting  the  oscillations  of  surface 
temperature  at  time  t  >■  0;  it  dies  away  as  t  increases,  leaving  the 
third  term  which  is  the  long  term  solution. 


A  truncated  cosine  function,  9l0),  can  also  be  used  to  represent 
temperature  variations  at  the  permafrost  table.  For  the  I'th  cycle. 


92if) 


|q  ip<t<ip+p\ 

AoCos{(a)-¥T^  lP^p\<t<lP-¥P2 
®  IP  +  /»i<f<(/+l)P 


(9) 


where  is  the  amplitude  cf  the  cosine  function  in  units  of  ®C,  Tm 
is  the  mean  value  of  the  cosine  function  during  the  entire  period  of 
P  and  is  not  in  general  equal  to  T.  Expanding  (9)  in  a  Fourier's 
series 


+  X  (OW-  «li) 

^  -l  (10) 


where  and  the  an  and  bn  are  the  Fourier 

coefficients  for  the  truncated  cosine  function. 
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(11) 


The  first  term  on  the  right  hand  side  of  (10)  is  T.  Then, 

T  ix,0)  =  T  +  GoX 


and 


T  (0,/)  =  7*  +  A„cos  iOM  -  O 

«-l  (12) 

The  steady-periodic  solution  for  (4),  (11)  and  (12)  can  be 

written  in  the  form 


T  (x,/)  =  r  +  C/<,X  +  X  ^H^'^COS  (CM  -  £«). 

n-1  (13) 


RISULTS  AKD  DISCUSSION 

The  first  term  on  the  right  hand  side  of  (8)  and  (13)  is  T.  The 
second  term  is  the  initial  temperature  variation  with  depth.  The 
third  term  is  the  temperature  oscillation  around  the  mean 
temperature  as  a  function  of  depth  and  time.  Figure  1  shows  the 
results  calculated  by  (8)  and  (13),  respectively  for  late  June.  The 
two  solutions  agree  at  depth  (greater  than  10  m  or  so)  but  generally 
disagree  near  the  permafrost  surface.  A  truncated  cosine  function 
is  a  better  fit  to  measured  temperatures  (Zhang,  1989)  so  that  it 
is  preferred  over  a  step  function  although  a  step  function  is 
simpler  to  use.  It  is  also  preferred  over  a  full  cosine  function 
which  introduces  too  much  warmth  during  the  summer  thaw  period. 

Figures  2  and  3  show  the  calculated  annual  oscillation  of 
temperature  at  various  depths  as  a  result  of  the  step  and  truncated 
cosine  boundary  conditions,  respectively.  The  step  and  truncated 
cosine  curves  in  Figures  2  and  3  are  the  boundary  conditions  with  n 
-  200  and  the  rest  of  the  curves  are  the  temperature  oscillation 
with  time  at  depths  from  2  to  10  m  with  an  increment  of  2  m  for  n  - 
100.  The  "wiggles"  in  the  step  curve  near  Pj  and  P2  in  Figure  2  or 

•  t 

and  ^2  in  Figure  3  are  caused  by  talcing  a  finite  number  of  terms 
in  the  series  in  (6)  and  (12)  . 

The  amplitude  of  the  temperature  oscillation  at  the  permafrost 
table  in  (8)  is 


Af  I 


(14) 
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and  in  (12)  it  is 


n«l 


(15) 


TEMPERftTUHE  fC) 
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Figure  1:  Temperature  variation  of  permafrost  with  depth  for  late 
June  calculated  for  solutions  with  a  step  (a)  and 
truncated  cosine  (b)  boundary  conditions  with 

Pi«70  days,  days,  P|*46  days,  days,  and  D-42.0  m^yr* 

starting  from  September  12. 


X  c*WnaF2ZT 

and  these  amplitudes  decrease  with  depth  according  to  m-i 
Since  D  and  <0  are  assumed  constant,  when  n  increases,  the  wave 

number  K  increases  and  wave  length  X»2mIK  decreases  until  the 
higher  harmonics  die  out  for  n  sufficiently  large.  The  higher 
harmonics  in  (8)  and  (13)  disappear  rapidly  moving  downwards  into 
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the  permafrost  and  the  step  and  truncated  cosine  functions  soon 
become  sinusoidal.  Figures  2  and  3  show  that  the  amplitude  of  the 
temperature  oscillation  decreases  with  depth  and  that  the  shape  of 
the  curves  approach  a  sinusoidal  function  at  depth.  Comparing 
Figure  2  with  Figure  3,  the  higher  harmonics  in  (13)  disappear  much 
faster  than  those  in  (8)  . 


Figure  2:  Annual  oscillation  of  temperature  at  2,  4,  6,  8,  and  10  m 
caused  by  a  step  boundary  condition  at  the  permafrost 
table  (x»0) . 


Figure  3:  Annual  oscillation  of  ten?)erature  at  2,  4,  6,  8,  and  10  m 
by  a  truncated  cosine  boundary  condition  at  the 
permafrost  table  (x»0) . 
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CONCLUSIONS 


The  thermal  regime  of  permafrost  within  the  depth  of  the  annual 
temperature  variations  is  often  analyzed  as  a  pure  heat  conduction 
problem  in  a  homogeneous  semi-infinite  region  with  the  boundary 
condition  a  sine  or  cosine  function  of  time  at  the  ground  surface  to 
represent  the  change  of  seasons.  However,  for  a  wet  active  layer, 
the  thermal  regime  of  permafrost  cannot  be  exactly  described  by  a 
sinusoidal  boundary  condition  at  the  ground  surface  due  to  phase 
change  in  the  active  layer  over  permafrost.  The  choice  of  the 
permafrost  table  as  the  upper  boundary  for  investigating  the  thermal 
regime  of  permafrost  overcomes  this  problem.  Temperature 
measurements  at  the  permafrost  table  show  that  the  temperature 
vari-ations  there  are  not  a  sinusoidal  function  of  time  since  the 
temperature  is  constrained  at  or  near  O^C  for  part  of  the  year.  The 
measured  data  show  that  the  temperature  variation  at  the  permafrost 
table  can  be  better  approximated  as  a  step  or  a  truncated  cosine 
function.  Analytical  solutions  to  the  heat  conduction  problem  were 
obtained  by  expanding  the  assumed  permafrost  surface  temperature 
variation  in  a  Fourier  series.  The  truncated  cosine  function  at  the 
surface  is  more  realistic  than  the  step  function.  When  n,  the 
number  of  terms,  is  large,  the  higher  harmonics  in  the  solutions  are 
less  important  and  disappear  rapidly  with  depth  in  the  permafrost. 
Both  the  step  and  truncated  cosine  function  at  the  permafrost 
surface  approach  a  sinusoidal  function  with  depth. 
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A  two-dimensional  finite  element  model  was  used  to  Investigate  the  thermal 
response  of  subsea  permafrost  and  gas  hydrates  to  changes  In  sea  level  and  climate 
over  a  121  Kyr  time  period  along  a  line  offshore  from  Lonely,  Alaska.  For  subsea 
permafrost  containing  brines,  the  spatial  distribution  of  the  Ice-bearlng 
permafrost  (IBP)  Is  predicted  to  lae  wedge  shaped  and  to  extend  only  19  km  offshore. 
There  Is  significant  lateral  heat  flow  throughout  the  IBP  section  offshore  and  the 
dept'a  to  the  IBP  table  Increases  almost  linearly  wltli  distance  offshore.  For 
subsea  permafrost  with  constant  thermal  parameters,  IBP  Is  predicted  to  extend  52 
km  offshore  (water  depth  -  50  m)  and  Is  nearly  Isothermal  beyond  4  km  offshore. 
Depth  to  the  IBP  table  Increases  almost  linearly  with  distance  offshore  and  then 
becomes  relatively  shallow  and  nearly  constant  In  dopth.  Seabed  temperatures  and 
the  assumed  sea  level  history  curve  are  especially  Important  In  determining  the 
current  distribution  of  subsea  permafrost.  The  full  thickness  of  IBP  onshore  can 
be  modeled  better  using  constant  thermal  parameters.  Depth  to  the  ILP  table 
offshore  may  be  modeled  better  when  It  Is  assumed  that  the  permafrost  contains 
brines.  The  predicted  depth  zone  for  stability  of  methane  gas  hydrates  Is  between 
220  m  and  650  m  near  shore.  For  subsea  perstafrost  containing  brines,  this  zone 
extends  32  km  offshore  compared  to  54  km  (to  the  55  m  water  depth)  for  constant 
thermal  parameters.  The  time  scale  for  producing  methane  gas  from  destabilized  gas 
hydrates  In  the  continental  shelf  near  Lonely  Is  on  the  order  of  10*  years,  much 
longer  than  previously  predicted. 


INTRODUCTION 

During  the  past  million  years,  there  Is 
evidence  for  repeated  glaciations  alx)ut  every  10^ 
years  (Shackleton  and  Opdyke,  1977).  The 
continental  shelves  of  the  Arctic  ocean  were 
emergent  during  these  glaciations,  up  to  current 
water  depths  of  120  m  (Bard  et  al.,  1990).  Cold 
climates  during  periods  of  emergence  favored  the 
formation  of  permafrost  and  the  stabilization  of 
gas  hydrates  In  the  shelves.  Rising  sea  levels 
during  Interglacial  periods  replaced  the  cold  air 
temperatures  on  the  shelves  by  much  warmer  sea 
water  temperatures.  As  a  result,  the  permafrost 
and  thawed  sediments  would  have  warmed  at  all 
depths,  and  permafrost  would  have  started  to  thaw 
from  )3oth  the  top  and  the  bottom.  Eventually,  gas 
hydrates  would  have  been  destabilized,  providing  a 
potential  source  of  methane  gas  to  the  atmosphere. 

Since  permafrost  thaws  very  slowly  In  response 
to  the  new  surface  boundary  conditions  after 
su)3mergence,  considerable  time  (on  the  order  of 
tens  of  millennia)  may  be  reciulred  to  completely 
thaw  the  permafrost.  Consequently,  Ice-bearlng 
permafrost  (IBP)  and  gas  hydrates  may  still  exist 
In  parts  of  the  continental  shelves  of  the  Arctic 
Ocean.  The  existence  of  subsea  permafrost  and  the 
presence  of  Ice  In  these  continental  shelves  has 
)7een  confirmed  by  drilling  In  shallow  water 
(usually  a  few  tens  of  meters  or  less  In  depth) 
off  the  coasts  of  the  USSR  (Molochushkln,  1978), 
Alaska  (Lewellen,  1973),  and  Canada  (Kackay, 

1972).  Subsea  permafrost  Is  also  Indicated  In  the 
Interpretations  of  offshore  Ismlc  data  (Hunter 
et  al.,  1976)  and  of  geophysical  logs  In  offshore 
petroleum  exploration  wells  (Osterkamp  et  al., 
1985) .  Destabilization  of  gas  hydrates  (by 
warming  the  sediments  In  the  continental  shelves) 
during  periods  of  high  sea  level  may  be  a  i^rlodlc 
source  of  atmospheric  methane  over  geological  time 
(Clarke  et  al.,  1986;  MacDonald,  1990).  Warming 
of  the  permafrost  and  sediments  and  permafrost 


thawing  eventually  causes  gas  hydrates  to  become 
unstaole  resulting  In  the  IDseratlon  of  large 
volumes  of  gas.  However,  the  permafrost  may  act 
as  a  seal  (t>ecause  of  Ice  In  the  pores)  preventing 
the  gases  from  escaping  until  sufficient  Ice  has 
been  thawed  to  generate  escape  routes  for  the  gas. 
MacDonald  (1990)  has  Investigated  the  time  scales 
for  the  response  of  the  permafrost  to  sutwergence 
using  a  one-dlmenslonal  analytical  model.  Since 
the  model  did  not  Include  latent  heat  effects,  the 
predicted  time  scales  are  expected  to  )3e  much  too 
short.  This  means  that  his  predictions  regarding 
the  time  scales  for  production  of  atmospheric 
methane  gas  by  destabilization  of  gas  hydrates  In 
continental  shelves  affected  )3y  permafrost  are  not 
correct . 

Calculations  using  one-dlmcnslonal  analytical 
B»dels  (e.g.  Mackay,  1972;  Lachenbruch  et  al., 

1982)  and  numerical  models  (Outcalt,  1985;  Nixon, 
1986)  for  the  transient  response  of  permafrost  to 
subeiergence  generally  suggest  the  presence  of 
relatively  thick  subsea  permafrost  even  In  areas 
of  deeper  water.  The  occurrence  of  subsea 
permafrost  Implies  the  potential  presence  of 
stability  zones  for  gas  hydrates.  However,  one- 
dlswnslonal  models  are  not  conpletely  satisfactory 
because  of  possible  lateral  heat  flow  In  the 
subsea  permafrost  particularly  near  shore.  In 
some  cases.  It  may  also  be  necessary  to  include 
distributed  latent  heat  effects  and  variable 
thermal  parameters  associated  with  the  potential 
presence  of  saline  pore  fluids  In  the  permafrost, 
which  do  not  appear  to  have  been  done. 

This  paper  presents  the  results  of  two- 
dimensional  numerical  modeling  of  the  thermal 
response  of  perstafrost  to  changes  In  sea  level  and 
climate,  which  Includes  the  effects  of  saline  pore 
fluids  on  latent  heat  and  thermal  parameters.  The 
region  considered  was  the  continental  shelf  of 
Alaska  near  Lonely.  Results  of  these  simulations 
are  used  to  evaluate,  at  the  present  time,  the 
spatial  dlstrllxitlon  of  IBP  In  the  continental 
shelf  near  Lonely  and  the  stability  zones  of  gas 


hydrates  chat  may  exist  within  or  under  th( 
permafrost . 

Additional  information  on  the  numerical 
simulations,  study  site,  choice  of  parameters, 
initial  conditions,  boundary  conditions,  and  other 
simulations  may  be  found  in  Fei  (1992). 

STUDY  SITS 

The  study  site,  offshore  from  Lonely  which  is 
about  135  Ion  southeast  of  Barrow,  Alaska,  was 
chosen  because  of  the  availability  of  data  from 
other  research  and  from  onshore  and  offshore 
petrole'uiTi  exploiacion  wells.  Available  data 
consist  of  geophysical  logs  and  samples  from  the 
J.W.  Dalton-1  (JWD)  well  onshore  and  the  Antares 
well  offshore  (Collett  et  al.,  1989;  Deming  et 
al.,  1992),  results  of  thermal  studies  in  five 
shallow  drill  holes  along  an  offshore  line  to  the 
northwest  from  Lonely  (Harrison  and  Osterkamp, 
1981),  and  results  of  bottom  sampling  and 
shoreline  erosion  studies  (Hopkins  and  Hartz, 

1978)  . 

The  onshore  surflcial  deposits  were  mapped  as 
interglacial  nearshore  and  lagoon  sand,  silty  fine 
sand,  and  pebbly  sand  (Hopkins  and  Hartz,  1978). 
Shallow  offshore  drilling  data  (Harrison  and 
Osterkamp,  1981)  Indicate  that  the  seabed 
sediments  are  fine-grained  to  about  the  30  m 
depth.  At  the  JWD  well,  the  deeper  well  logs 
indicate  relatively  coarse  material  (conglomerate 
and  sandstone)  down  to  the  270  m  to  300  m  depth  or 
so  overlying  finer  material  (siltstone)  (Collett 
et  al.,  1989).  There  is  a  change  in  the 
temperature  gradient  in  this  depth  interval 
(Lachenbruch  and  Marshall,  1986).  We  interpret 
the  base  of  the  IBP  to  be  at  360  m  to  366  m  (about 
-2.2°C)  where  slight  changes  in  the  resistivity 
log  and  the  temperature  gradient  occur.  Other 
onshore  wells  in  this  region  Indicate  a  similar 
lithology  but  generally  with  somewhat  finer 
material . 

The  nearest  offshore  well  (Antares,  EXXON)  is 
at>out  24  km  distant  on  a  line  laearing  N55*E  from 
Lonely  in  alxjut  15  m  of  water.  Well  logs 
indicate  that  relatively  fine-grained  material 
exists  in  the  upper  section  (above  190  m)  where 
permafrost  might  be  expected  with  some  coarser 
material  from  190  m  to  312  m.  The  use  of 
geophysical  logs  to  determine  the  presence  or 
absence  of  IBP  is  very  difficult  because  of  the 
lack  of  contrast  in  physical  properties  between 
the  thawed  material  and  rny  warm  and  marginally 
ice-bearing  permafrost  which  would  be  thawing  from 
)xith  the  top  and  (xittom.  According  to  the  gamma 
ray  log,  there  are  two  similar  sections  Just  a)30ve 
272  m  and  305  m.  Increases  in  the  resistivity  log 
(DID  above  272  m  but  not  akxjve  305  m  suggest  that 
IBP  may  exist  a)Dove  272  m  which  is  our 
interpretation  of  the  log.  However,  other  factors 
(e.g.,  changes  in  pore  fluid  salinity)  besides 
perm.afrost  could  produce  these  increases. 

The  harmonic  mean  thermal  conductivity 
determined  from  drill  cuttings  from  the  JWD  well 
in  the  frozen  interval  from  113  m  to  269  m  is  2.57 
W(mK)'i  (Deming  et  al.,  1992).  Values  of  2.5 
W(mK)'^  and  1.7  W(mK)''  were  used  for  the  frozen 

and  unfrozen  material  respectively.  The  porosity 
was  assumed  to  be  0.2.  For  fine-grained 
sediments  containing  brines,  the  temperature- 
dependent  thermal  parameters  were  calculated  using 
an  unfrozen  water  content,  0u  =  AT®,  where  T  is 
the  magnitude  of  the  temperature  (°C)  and  A  and  B 
are  empirical  constants  (A  =  0.1435,  B  =  -0.902) 
(Fei,  1992).  Freezing  point  depression  was 
assumed  to  be  -1.63“C.  The  volumetric  heat 
capacity  and  the  thermal  diffuslvity  were 
calculated  according  to  the  methods  used  by 
Osterkamp  (1987). 


Deming  et  al.  (1992)  suggest  a  large  heat  flow 
of  ■  80  mW-m"2  for  the  JWD  well.  However,  the 
average  measured  temperature  gradient  and  the 
above  frozen  conductivity  yield  a  more  normal 
value  of  "  54  mW-m'^  in  the  upper  interval.  For 
the  simulations,  a  compromise  value  of  65  mW-m'7 
was  used. 

Information  on  the  near  surface  seabed 
sediments,  temperatures,  and  depth  to  IBP  was 
obtained  in  five  shallow  drill  hole.s  which  were 
rotary  jet  drilled  along  a  line  bearing  N32°w  near 
the  DEW  site  at  Lonely  (Harrison  and  Osterkamp, 
1981)  .  Distances  from  shore  were  paced  and  may 
contain  significant  errors. 

Figure  1  shows  temperature  profiles  measured 
in  these  holes,  two  to  three  weeks  after  drilling, 
which  are  thought  to  be  within  a  few  hundredths 
degree  of  equilibrium  values.  The  profiles  show  a 
decreasing  thermal  gradient  with  distance  (time) 
offshore.  Assuming  the  current  average  shoreline 
erosion  rate  (4.7  ma'i),  the  hole  at  7770  m 
offrhore  has  been  submerged  for  a)Dout  seventeen 
centuries.  Approximate  mean  seabed  temperatures 
can  be  obtained  by  projecting  the  thermal 
gradients  in  the  deeper  part  of  the  holes  up  to 
the  seabed.  However,  the  presence  of  a  phase 
(boundary  often  makes  it  difficult  to  do  this 
reliably. 


Figure  1.  Subsea  permafrost  temperatures  measured 
in  five,  shallow,  offshore  holes  along  a  line 
bearing  N32'W  at  Lonely,  Alaska.  The  numbers  next 
to  the  profiles  are  the  estimated  distances 
offshore  which  were  paced  and  may  contain 
significant  errors. 


In  warm  subsea  permafrost  which  contains 
brines,  the  amount  of  ice  appears  to  increase 
gradually  with  depth  (Osterkamp  and  Hat  risen, 

1982;  Osterkamp  et  al.,  1989).  Consequently, 
there  is  no  distinct  change  in  temperature 
gradient  nor  properties  at  the  IBP  table  but  only 
a  gradual  curvature  of  the  temperature  profile 
which  makes  it  difficult  to  detect  the  presence  of 
ice  from  the  drllllrg  data  or  the  temperature  data 
(Fig.  1).  Therefore,  the  holes  were  heated  to 
determine  the  presence  of  ice  using  the  method 
described  )3y  Osterkamp  and  Harrison  (  1980),  Our 
interpretation  of  the  results  of  drilling  and 
heating  the  holes  indicate  that  the  IBP  table 
occurs  between  6  m  and  15  m  below  the  seabed  along 
this  line  (Table  1). 


Table  1,  Data  for  holes  drilled  at  Lonely  along  a  line  N32°W  starting  at  a  point  on  the  shoreline  which 
was  N66°W  from  the  DEW  site  radar  dome.  The  hole  number  is  the  distance  from  shore  in  meters  determined 
by  pacing. 
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7,70  m 


Simulations  were  done  for  two  lines,  one 
bearing  N32°W  (Table  1)  which  crosses  the  shelf  at 
an  angle  and  one  bearing  N20°E  which  is 
perpendicular  to  the  depth  contours  and  passes 
about  15  )tm  to  the  west  of  the  Antares  well. 

Seabed  profiles  for  these  lines  were  constructed 
using  results  from  the  above  studies  and 
bathymetric  maps.  Two  simulations  for  the  latter 
line  are  reported  herein,  one  for  sediments  with 
fresh-water  pore  fluids  and  one  for  sediments 
containing  sea  water  brines.  Our  strategy  was  to 
perform  two  simulations  for  two  near-extreme  cases 
with  regard  to  their  effects  on  the  thermal  state 
of  the  permafrost.  Additional  details  and 
simulations  are  discussed  in  Fei  (1992). 


The  initial  thermal  conditions  are  unicnown. 
However,  the  simulations  were  started  using  steady 
state  conditions  (Lachenbruch,  1957)  far  enough 
into  the  past  (121  Kyr  BP)  to  allow  any  transients 
(Oster)camp  and  Gosin)c,  1991)  associated  with  the 
initial  conditions  to  disappear.  The  surface 
boundary  condition  depends  on  whether  the  surface 
is  emergent  or  submergent.  For  the  emergent 
surface  temperature,  the  paleocllmatlc  temperature 
history  of  Maximova  and  Romanovs)ry  (1988), 
modified  for  the  Lonely  region  (Oster)camp  and 
Gosin)t  1991),  was  used.  Their  history  gives  good 
agreement  with  the  current  permafrost  thic)tnesa  at 
Prudhoe  Bay  (Oster)camp  and  Goaln)c,  1991).  Since 
there  does  not  appear  to  be  a  sea  level  history 
curve  for  this  portion  of  the  Beaufort  Sea  shelf, 
the  timing  tor  emergence  and  submergence  was 
obtained  from  the  sea  level  history  curve  of  Bard 
et  ai.  (  1990)  modified  to  ta)ce  shoreline  erosion 
into  account  (Fei,  1992).  This  procedure  produces 
depths  to  the  seabed  that  are  somewhat  shallower 
than  observed  up  to  20  )an  offshore.  For  the 
submerged  surface,  current  seabed  temperatures 
determined  from  measurements  in  shallow  drill 
holes  at  Prudhoe  Bay,  Barrow,  and  Lonely,  adjusted 
for  water  depth,  were  used  (Fig.  1  and  Osterlcamp 
and  Harrison,  1982,  1985).  Additional  information 
is  provided  by  Fei  (1992). 


Figure  2  shows  the  predicted  current 
temperature  distribution  in  the  continental  shelf 
along  the  line  bearing  N20*E  from  Lonely.  The 
uneven  character  of  the  isotherms  in  Figure  2  is 
caused  by  the  plotting  algorithm.  IBP  (-1.6*C 
isotherm)  is  predicted  to  extend  only  about  19  )an 
offshore.  This  surprising  prediction  is  a  result 
of  the  relatively  large  heat  flow  and  low  ice 
content  assumed  for  the  IBP. 


D4«aBet(km) 

Figure  2.  Predicted  offshore  temperature 
distribution  and  stability  zone  for  gas  hydrates 
at  the  present  time  near  Lonely  for  fine-grained 
sediments  containing  brines.  The  ice-bearing 
permafrost  is  defined  by  the  -1.6°C  isotherm. 


Lateral  heat  flow  is  large  in  the  IBP  near  the 
coast,  but  it  is  also  significant  in  the  rest  of 
the  IBP  and  in  the  thawed  material  under  it  and 
seaward  of  it.  Beyond  50  )ui\  offshore,  at  deeper 
depths  the  curved  isotherms  and  associated  lateral 
heat  flow  appear  to  be  the  result  of  recently 
thawed  IBP  near  the  seabed.  These  results  suggest 
that  one-dimensional  thermal  models  that  assume 
there  is  no  lateral  heat  flow  may  not  correctly 
predict  the  thermal  regime  of  subsea  permafrost 
containing  brines. 

The  predicted  depth  to  the  base  of  IBP  onshore 
is  318  m  compared  to  the  observed  depth  of  366  m, 
a  difference  of  13%.  In  the  offshore  region,  the 
base  of  the  IBP  rises  rapidly  with  distance 
offshore.  A  comparison  with  the  Antares  well  is 
difficult  since  it  is  about  15  )tm  east  of  this 
profile  in  a  water  depth  of  15  m.  If  water  depth 
is  used  as  a  criterion,  this  calculated  profile 
predicts  an  absence  of  IBP  at  the  Antares  well. 
However,  if  distance  offshore  is  used,  IBP  at  the 
Antares  well  would  be  predicted  to  a  depth  cf 
about  200  m  compared  to  the  272  m  observed. 

Depth  to  the  IBP  table  increases  almost 
linearly  with  distance  offshore  to  about  66  m 
below  the  seabed  at  19  itm  offshore.  At  7.8  )cm 
offshore,  the  predicted  thic)cness  of  the  thawed 
layer  at  the  seabed  is  21  m.  At  the  same  distance 
offshore,  along  the  line  bearing  N32°W,  the 
observed  thic)tness  is  7  m  to  15  m  (Table  1).  I: 
water  depth  were  used  as  a  criterion,  the 
difference  between  observed  and  calculated  values 
would  be  greater. 


Fine-grained  Material  with  Constant  Thermal 
Properties 

Figure  3  shows  Che  predicted  current 
temperature  distribution  along  Che  line  bearing 
N20°E  from  Lonely.  A  nearly  isothermal  IBP  section 
extends  from  4  km  to  52  km  offshore  where  Che 
water  depth  is  about  50  m.  Up  Co  about  24  km 
offshore,  the  IBP  section  is  wedge-shaped  and 
beyond  24  km  it  is  nearly  tabular  with  a 
relatively  flat  vertical  tip.  The  relatively 
thick  tabular  section  appears  Co  be  a  result  of 
Che  assumed  sea  level  curve,  This  curve  indicates 
chat  depths  of  less  Chan  50  m  were  continuously 
exposed  to  cold  air  temperatures  for  more  than  70 
kyr  and  only  submerged  during  Che  last  11  kyr. 

This  would  produce  a  chick  section  of  IBP  out  to 
the  current  50  m  isobath. 
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Figure  3.  Predicted  offshore  temperature 
distribution  and  stability  zone  for  gas  hydrates 
at  the  present  time  near  Lonely  for  fine-grained 
sediments  with  constant  thermal  parameters.  The 
ice-bearing  permafrost  is  defined  by  the  -1.6*C 
isotherm . 


Lateral  heat  flow  is  large  within  4  tan  of 
shore  and  in  the  thawed  material  beyond  and  under 
Che  tip  of  Che  IBP.  Predicted  depth  to  the  base 
of  Che  IBP  onshore  (355  m)  is  close  to  the 
observed  value  (3£6  mj .  This  suggests  that  the 
full  thickness  of  the  IBP  iriay  be  modeled  better 
with  constant  thermal  properties  than  with 
temperature-dependent  properties  and  distributed 
latent  heat  effects  associated  with  the  presence 
of  brines  in  the  permafrost.  A  similar  result  was 
found  by  Oscerkamp  and  Goslnk  (1991)  when  modciing 
changes  in  permafrost  thickness  at  Prudhoe  Bay  in 
response  to  changes  in  paleocl imate .  In  the 
offshore  region,  the  predicted  base  of  the  IBP 
rises  rapidly  to  24  km  offshore  and  then  more 
slowly  farther  offshore.  Comparison  with  the 
Ancares  well  is  again  difficult  .  At  the 
ecpuivalenc  water  depth  (15  m) ,  the  predicted  base 
of  die  IBP  would  be  at  about  210  m  below  sea 
level,  and  at  the  equivalent  distance  offshore, 
about  270  m  compared  to  Che  observed  depth  of  272 
m. 

Depth  CO  the  ISP  table  increases  almost 
linearly  with  distance  offshore  and  reaches  a 
maximum  value  of  66  m  below  the  seabed  at  22  km 
offshore.  It  Chen  rises  to  within  30  m  of  the 
seabed  and  remains  at  nearly  a  constant  depth  to 
the  tip  of  Che  IBP.  This  behavior  Is  the  combined 
result  of  the  sea  level  history  curve  and  seabed 
temperatures.  Sea  levels  rose  rapidly  from  the 
glacial  minimum  less  chan  20  kyr  BP  to  about  4  kyr 


BP.  As  a  result,  Che  shelf  was  rapidly  covered 
with  deep  cold  water  (about  -1.5°C),  In  the  last 
4  kyr,  submergence  (due  Co  sea  level  rise  and 
shoreline  erosion)  occurred  in  shallow  warm  water 
(alDouc  -1°C)  .  In  Che  shallow  water,  Che  larger 
temperature  gradient  in  Che  chawed  material 
between  Che  seabed  and  Che  IBP  Cable  (-1.63°C) 
allows  for  greater  heat  flow  and  therefore  a 
faster  thawing  race.  The  predicted  thickness  of 
Che  chawed  layer  below  the  seabed  is  38m  at  7.8  km, 
offshore,  much  greater  chan  observed  along  our 
line  of  shallow  holes  (Table  1.).  While  Chawing 
near  the  se?bed  for  this  case  is  initially  greater 
chan  when  brines  are  present,  by  19  km  offshore 
Che  thicknesses  of  the  thawed  layers  are  about 
equal . 

Simulations  for  looth  materials  were  carried 
out  along  our  line  of  shallow  holes  (N32°W)  using 
Che  aloove  parameters  and  conditions  except  for  a 
larger  geothermal  heat  flow  of  80  mW-m"^  at 
suggested  by  Deming  ec  al.  (1992) .  In  these 
cases,  predictions  for  the  depth  to  the  base  of 
IBP  onshore  were  much  coo  shallow  suggesting  that 
this  choice  of  heat  flow  may  be  coo  large  for 
Lonely . 

Stability  Zone  for  Gas  Hydrates 

Gas  hydrates  are  stable  over  a  limited  range 
of  pressures  and  temperatures  chat  can  be  found  in 
association  with  permafrost,  including  subsea 
permafrost  (e.g.  Kvenvolden  and  McMenamin,  1980: 
Collect  ec  al.,  1988).  The  stable  region  can  be 
determined  from  the  phase  equilibrium  diagram 
(Sloan,  1990) .  It  was  assumed  chat  Che  only  gas 
in  the  hydrate  was  methane  and  chat  pressures 
could  be  converted  to  depth  using  hydrostatic 
conditions  in  the  sediments.  Latent  heat  effects 
due  CO  formation  or  decay  of  the  gas  hydrates  were 
neglected,  with  these  assumptions,  the  current 
approximate  stability  regions  for  methane  gas 
hydrates  in  the  continental  shelf  near  Lonely  were 
mapped  on  the  two-dimensional  depth-temperature 
diagrams  in  Fig.  2  and  3. 

In  fine-grained  sediments  containing  brines, 
gas  hydrates  can  potentially  exist  up  to  32  km 
offshore.  The  stability  zone  is  below  220  m  and 
above  650  m  near  shore,  about  the  same  as  Chat 
shown  in  Lachenbruch  ec  al.  (1988)  for  the  JWD 
well.  For  constant  thermal  parameters,  the  depth 
range  of  the  stability  zone  near  shore  is  similar: 
however,  the  stability  zone  extends  up  to  about  54 
km  offshore  where  the  water  depth  is  about  55  m. 
Using  the  sea  level  history  curve  of  Bard  ec  al. 
(1990),  depths  of  50  m  to  55  m  were  submerged 
aloout  10*  years  ago.  Thus,  the  time  scale  for 
producing  methane  gas  from  doscablllzed  gas 
hydrates  in  the  continental  shelf  near  Lonely  (for 
constant  thermal  parameters)  is  about  four  times 
greater  than  predicted  by  MacDonald  (1990). 

summary: 

A  two-dimensional  finite  element  model  was 
used  to  evaluate  the  thermal  response  and  current 
thermal  regime  of  permafrost  to  changes  in  sea 
level  and  climate.  This  information  was  used  to 
calculate  the  current  stability  zone  for  gas 
hydrates  in  the  continental  shelf.  ■  study 
site,  offshore  from  Lonely,  was  chosv,  because  of 
the  availability  of  data  from  ocher  research  and 
from  two  oil  exploration  wells,  one  onshore  and 
one  offshore.  These  data  were  used  to  provide 
information  on  sediment  types,  depths  to  the 
permafrost  cable  and  base,  thermal  properties  of 
the  sediments,  heat  flow,  shoreline  erosion  rate, 
and  boundary  conditions  for  modeling.  Simulations 
were  carried  out  for  two  extreme  cases,  permafrost 
with  constant  thermal  parameters  and  permafrost 
containing  sea  water  brines  which  introduce 
distributed  latent  heat  effects  and  cause  the 


thermal  parameters  to  have  a  strong  dependence  on 
temperature.  These  two  cases  span  a  wide  range  of 
thermal  parameters  and  produce  results  for  the 
thermal  states  of  the  permafrost  which  are  quite 
different . 

For  the  case  when  the  permafrost  contains 
brines,  ice-bearing  permafrost  (IBP)  is  predicted 
to  extend  only  19  l<m  offshore.  This  is  a  result 
of  the  relatively  large  heat  flow  (65  mW-m"2)  and 
low  ice  content  (porosity)  assumed  for  the  IBP. 
Lateral  heat  flow  is  significant  throughout  the 
IBP  and  in  the  thawed  material  under  it  and 
seaward  of  it.  This  suggests  that  one-dimensional 
thermal  models  may  not  correctly  predict  the 
thermal  regime  of  subsea  permafrost  containing 
brine.  Depth  to  the  IBP  table  increases  almost 
linearly  with  distance  offshore  and,  at  7.8  loti,  is 
28  m  compared  to  the  7  m  to  15  m  observed  along 
another  line  at  the  same  distance  offshore. 

For  the  case  when  thermal  properties  are 
constant,  IBP  is  predicted  to  extend  52  lun 
offshore  where  the  water  depth  is  50  m.  The  IBP 
is  primarily  isothermal,  wedge  shaped  to  24  km 
offshore,  nearly  tabular  beyond  24  km  offshore, 
and  has  a  relatively  flat  vertical  tip.  Depth  to 
the  IBP  table  increases  almost  linearly  with 
distance  offshore  over  the  wedge  shaped  section 
and  is  relatively  shallow  and  nearly  constant  over 
the  tabular  section.  This  behavior  is  the 
combined  result  of  the  assumed  sea  level  history 
curve  and  seabed  temperatures.  At  7.8  km 
offshore,  the  predicted  thickness  of  the  thawed 
layer  at  the  seabed  is  38  m  compared  to  7  m  to  15 
m  observed  along  another  line  at  the  same  distance 
offshore . 

Predicted  depth  to  the  base  of  IBP  onshore 
obtained  by  applying  the  model  with  constant 
thermal  properties  (355  m)  is  in  better  agreement 
with  the  observed  value  (366  m)  compared  to  the 
prediction  using  temperature-dependent  thermal 
properties  (318  m) .  This  suggests  that  Che  full 
thickness  of  the  IBP  can  be  modeled  better  with 
constant  thermal  properties.  However,  the  depth 
to  Che  IBP  table  appears  to  be  predicted  better 
using  the  model  where  the  permafrost  contains 
brines.  While  chawing  of  subsea  permafrost  is  a 
complex  process  possibly  involving  convective  heat 
flow  by  movement  of  the  brines,  this  suggests  chat 
most  of  the  deeper  permafrost  may  be  relatively 
free  of  brines  but  chat  the  near-surface 
permafrost  contains  brines,  although  ocher  factors 
may  play  a  role. 

Simulations  using  the  large  heat  flow  (80  mW- 
m"-)  suggested  by  Deming  ec  al.  (1992)  resulted  in 
predictions  for  the  base  of  the  IBP  which  were 
much  too  shallow. 

The  predicted  zone  for  stability  of  methane 
gas  hydrates  is  below  about  220  m  and  above  about 
650  m  near  shore.  When  the  permafrost  contains 
brines,  this  stability  zone  extends  32  km 
offshore.  For  permafrost  with  constant  thermal 
properties,  the  stability  zone  extends  54  km 
offshore  to  the  55  m  water  depth.  Submergence 
occurred  at  this  water  depth  about  10*  years  ago. 
This  indicates  that  the  time  scale  for  producing 
methane  gas  from  destabilized  gas  hydrates  in  the 
continental  shelf  near  Lonely  is  about  four  times 
greater  chan  predicted  by  MacDonald  (1990). 
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Modeling  the  Response  of  Permafrost  and  Gas  Hydrates  to  Changes  in  Sea  Level 
and  Climate 

T.  E.  Qsterkamp  and  T.  Fei  (Both  at:  Geophysical  Institute,  University  of  Alaska 

Fairbanks,  Fairbanks,  AK  99775-0800;  907-474-7548) 

A  two-dimensional  finite  element  model  was  used  to  investigate  the  thermal  response 
of  subsea  permafrost  and  gas  hydrates  to  changes  in  sea  level  and  climate  along  a  line 
offshore  from  Lonely,  Alaska.  Simulations  were  carried  out  for  two  extreme  cases: 
permafrost  with  constant  thermal  parameters;  and,  permafrost  with  temperature-depen¬ 
dent  thermal  parameteis  and  distributed  latent  heat  effects  due  to  the  presence  of  bi'ines. 
For  subsea  permafrost  containing  brines,  the  ice-bearing  permafrost  (IBP)  is  predicted  to 
be  wedge  shaped  and  to  extend  only  19  km  offshore.  There  is  significant  lateral  heat  flow 
throughout  the  IBP  section  and  the  depth  to  the  IBP  table  increases  almost  linearly  with 
distance  offshore.  For  subsea  permafrost  with  constant  thermal  parameters,  IBP  is  pre¬ 
dicted  to  extend  52  km  offshore  (water  depth  -  50  m)  and  is  nearly  isothermal  beyond 
4  km  offshore.  It  is  wedge  shaped  up  to  24  km  offshore,  tabular  beyond,  and  has  a  flat 
vertical  tip.  Depth  to  the  IBP  increases  almost  linearly  with  distance  offshore  over  the 
wedge-shaped  section  and  is  relatively  shallow  and  nearly  constant  over  the  tabular 
section.  Seabed  temperatures  and  the  assumed  sea  level  history  curve  are  especially 
important  in  determining  the  current  distribution  of  subsea  permafrost  The  full  thickness 
of  IBP  onshore  can  be  modeled  better  with  constant  thermal  parameters  compared  to 
permafrost  that  contains  brines.  Depth  to  the  IBP  table  offshore  may  be  modeled  better 
when  it  is  assumed  that  the  permafrost  contains  brines.  The  predicted  zone  for  stability  of 
methane  gas  hydrates  is  between  220  m  and  650  m  near  shore.  For  subsea  permafrost 
containing  brines,  this  zone  extends  32  km  offshore  compared  to  54  km  (to  the  55  m 
water  depth)  for  constant  thermal  parameters.  The  time  scale  for  producing  methane  gas 
from  destabilized  gas  hydrates  in  the  continental  shelf  near  Lonely  is  on  the  order  of  10* 
years,  much  longer  than  previously  predicted. 
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Modeling  the  Response  of  Permafrost  to  Surface  Temperature  Changes  in  the 
Alaskan  Arctic 

T.  Zhang  and  T.  E.  Osterkamp  (Both  at;  Geophysical  Institute,  University  of  Alaska 

Fairbanks.  Fairbanks,  AK  99775-0800;  907-474-7548) 

Climatological  data  and  modeling  have  been  used  to  investigate  the  response  of  perma¬ 
frost  in  the  Alaskan  Arctic  north  of  the  Brooks  Range  to  changes  in  climate  over  the  last 
century  or  so.  Air  temperatures  in  this  region  are  strongly  correlated,  vary  with  periods  of 
10  years  and  47  years,  and  follow  the  same  general  trends  as  the  North  American  Arctic. 
Precipitation,  primarily  snowfall,  which  was  heavier  during  cold  periods  and  lighter 
during  warm  periods,  was  poorly  correlated  between  stations,  and  showed  a  periodicity 
of  37  years  at  Barrow  and  34  years  at  Barter  Island.  These  results  suggest  that  the  effects 
of  changes  in  air  temperatures  on  permafrost  temperatures  may  have  been  modified  by 
changes  in  snow  cover.  A  numerical  model  using  North  American  Arctic  air  tempera¬ 
tures  at  the  permafrost  surface  showed  that  the  permafrost  temperature  variations  de¬ 
pended  strongly  on  the  choice  of  initial  conditions,  that  air  temperature  changes  were 
sufficiently  large  to  account  for  the  observations,  and  that  air  temperature  changes  pre¬ 
dated  permafrost  temperature  changes.  An  alternative  explanation  of  the  last  result  is  that 
the  shallow  permafrost  temperatures  were  colder  than  the  assumed  initial  conditions.  The 
model  with  Barrow  air  temperanires  (1923-1991)  applied  at  the  ground  surface  with  no 
snow  cover  predicted  a  slight  cooling  of  permafrost  temperatures  at  present  Adding  the 
snow  cover  caused  the  permafrost  to  warm  with  about  the  same  magnitude  and  penetra¬ 
tion  depth  as  observed.  The  results  may  explain  why  different  sites  show  a  warming, 
little  or  no  change,  or  a  cooling  of  permafrost  temperatures.  Current  interpretations  of 
observed  permafrost  temperature  changes  are  unsatisfactory.  The  data  indicate  that  the 
warming  at  most  sites  occurred  just  before  or  later  than  1920.  However,  air  temperature 
trends  from  about  1920  to  present  appear  to  have  been  within  a  degree  or  so  of  the  long¬ 
term  mean  at  Barrow  (-12.5^0.  It  is  difficult  to  see  how  air  temperature  changes  alone 
could  have  proorced  the  observed  warming  at  Prudhoe  Bay  (2”C  to  3*C).  Current  differ¬ 
ences  between  air  and  permafrost  temperatures  are  about  3”C  to  6”C.  If  an  equilibrium 
initial  condition  is  assumed  for  the  permafrost  in  1923,  then  this  difference  must  have 
been  on  the  order  of  a  degree  or  less.  This  implies  a  very  thin  snow  cover  prior  to  that 
time,  non-equilibriui conditions  in  1880,  or  some  other  unknown  factor.  Variations  in 
snowfall  may  be  large. y  responsible  for  the  observed  wanning  of  permafrost  tempera¬ 
tures;  however,  there  is  insufficient  long-term  data  to  directly  test  this  hypothesis. 
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Considerations  in  Determining  Thermal  Diffusivity 

from  Temperature  Time  Series  Using  Numerical  Methods 

T.  Zhang  and  T.  E.  Osterkamp 

Geophysical  Institute,  University  of  Alaska  Fairbanks,  99775-  0800 

ABSTRACT 

This  paper  investigates  numerical  methods  for  determining  the  in  situ  appar¬ 
ent  thermal  diffusivity,  D,  of  the  active  layer  and  permafrost  from  repeated 
measurements  of  vert’cal  temperature  proSles.  The  results  can  be  applied 
to  any  soil  or  rock  where  heat  Sow  is  conductive.  An  analytical  expres¬ 
sion  was  derived  for  D  when  unfrozen  water  is  present  in  frozen  soils  and 
permafrost.  The  usual  numerical  expression  for  D  (termed  model  I)  was  ex¬ 
tended  to  include  higher  order  terms  (model  II).  This  substantially  reduces 
truncation  errors  but  increases  measurement  errors  somewhat.  Extremely 
accurate  temperature  measurements  (approaching  ±0.0VC)  are  required  to 
apply  the  method.  Synthetic  temperature  time  series  were  used  to  evaluate 
and  compare  predictions  of  models  I  and  II.  Model  I  produced  spikes  (large 
positive  and  negative  values)  in  the  calculated  D  at  times  where  the  ground 
temperatures  were  near  a  maximum  or  minimum.  Using  model  II,  the  spikes 
disappeared  or  were  substantially  reduced  in  magnitude  and  duration  and 
errors  in  D  at  other  times  were  reduced  compared  to  model  I.  Application  of 
the  method  requires  acquisition  of  very  precise  temperature  measurements  at 
appropriate  space  (Ax)  and  time  (At)  intervals  at  a  depth  where  measurable 
temperature  changes  occur.  Selection  of  values  for  Ax  and  At  must  take  into 
account  the  accuracy  of  the  temperature  measurements,  duration  and  ampli¬ 
tude  of  the  temperature  changes,  depth,  and  the  expected  values  for  D.  In 
general.  Ax  should  be  much  less  than  the  depth  of  interest  and  At  must  less 
than  the  period  or  duration  of  surface  temperature  changes.  For  the  active 
layer.  Ax  can  range  from  a  few  centimeters  to  a  tenth  of  a  meter  or  so  and  At 
from  a  few  minutes  to  an  hour  or  so.  For  the  near-surface  permafrost  within 
the  depth  of  annual  temperature  variations.  Ax  can  range  from  a  few  tenths 
meter  to  about  one  meter  and  At  from  one  day  to  sev^T^J  weelcs.  Below  the 
depth  of  annual  temperature  variations  in  the  permafrost,  values  for  Ax  of 
several  meters  and  At  of  a  year  or  more  may  be  appropriate. 
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INTRODUCTION 


Investigations  of  the  thermal  regime  of  the  active  layer  and  perma&ost  require 
an  understanding  of  their  thermzd  properties,  particularly  thermal  conductivity,  vol¬ 
umetric  heat  capacity  and  thermal  diifusivity.  The  thermal  diffusivity  is  of  primary 
importance  in  natural  systems  since  it  controls  the  rate  of  thermal  response.  McGaw 
et  al.,  (1978)  developed  a  method  to  determine  thermal  diffusivity  in  situ  for  the  active 
layer  and  perm2drost  from  a  temperature  time  series.  This  method  has  also  been  used 
by  Nelson  et  al.  (1985),  Outcalt  and  Hinkel  (1989,  1990),  Zhang  (1989),  and  Hinkel 
et  al.  (1990).  While  it  appears  that  good  results  can  be  obtained,  reported  in  situ  val¬ 
ues  for  thermal  diffusivity  are  sometimes  large,  contain  zero  values,  may  be  negative, 
and  have  considerable  scatter.  Using  a  synthetic  thermal  profile,  Hinkel  et  al.  (1990) 
showed  that  large  positive  and  negative  values  (“spikes”)  were  a  result  of  the  method 
and  occur  when  the  rate  of  change  of  the  temperature  gradient  becomes  small.  They 
attributed  other  large  fluctuations  in  thermal  diffusivity  to  non-conductive  heat  flow 
processes  in  the  active  layer;  a  conclusion  which  is  rendered  imcertain  by  the  variable 
results  for  D. 

This  paper  investigates  further  the  application  of  numerical  methods  for  deter¬ 
mining  the  thermal  diffusivity  of  the  active  layer  and  permafrost  &om  temperature 
time  series. 


REVIEW  OF  THEORY 

The  one-dimensional  transient  heat  conduction  equation  is 


dx 


(1) 


where  K  is  the  thermal  conductivity  (Wm“*®C~*),  T  is  the  temperature  (®C),  x  is 
the  distance  coordinate  (m),  Cv  is  the  volumetric  heat  capacity  (Jm~’*’C~‘)  and  t  is 
the  time  (sec). 

In  the  application  of  (1),  it  is  often  assumed  that  K  does  not  depend  on  position 
(i.e.,  the  soil  is  homogeneous).  In  a  thawed  active  layer  and  in  a  frozen  active  layer 
and  permafrost  that  does  not  contain  uxifrozen  water,  K  and  Cv  are  approximately 
independent  of  temperature.  With  these  assumptions,  (1)  becomes 


ax2  ~  D  dt 


(2) 
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where  D  =  K/Cv  is  the  thermal  difFusivity  {m^/sec). 

Equation  (2)  has  been  used  in  the  past  to  calculate  an  apparent  thermal  difFu¬ 
sivity  D  =  Tt/Txz,  where  Tt  and  T**  are  the  numerical  approximations  of  dTjdt  and 
d^T/dx^,  respectively.  This  relation  has  a  singular  point  at  Ttt  =  0  which  may  create 
difficulties  when  using  it  to  calculate  D. 

If  the  active  layer  and  permafrost  are  partieJly  frozen  zmd  contain  unfrozen  water, 
then  K  and  Cv  are  functions  of  temperature.  Equation  (1)  becomes  (Guymon,  1984) 

4.  =  r  ^ 

dz^  dT\dx)  "  dt  "p,dt 

where  L  is  the  volumetric  latent  heat  of  fusion  9i  is  the  volumetric  ice  content 

(m^/m^),  and  pi  and  pu  are  the  densities  of  ice  and  unfrozen  water,  respectively. 
Equation  (3)  can  be  written 


d^T  dK(dTV_  dT 

dT\dx) 


(4) 


where  Ca  =  Cv  +  L{dB^(dT)  is  the  apparent  volumetric  heat  capacity  and  is  the 
volumetric  unfrozen  water  content  (m’/m^)  (eg.  Osterkamp,  1987).  An  apparent 
thermal  diffu  ivity  for  frozen  soils  and  permafrost  containing  unfrozen  water  can  be 
defined  as  follows.  For  saturated  soils,  common  in  permafrost  and  the  active  layer. 


K  =  "  (5) 

where  K,,  Ki  and  Ku  ve  the  thermal  conductivities  of  the  soil,  ice  and  unfrozen  water 
components,  and  is  the  porosity  in  the  frozen  state.  The  temperature  dependence  of 
the  thermal  conductivity  of  the  constituents  is  ignored  over  the  range  of  temperatures 
usually  found  in  the  field.  Usings.  =  p^w^lpu  <snd  u*.  =  A(-T)^,  where  p^  is  the  bulk 
density  of  the  soil  {kg/m^),  p^  is  the  density  of  unfrozen  water,  is  the  gravimetric 
water  content  {kg/kg),  and  A  and  B  are  constants,  then 

If  =  -Ri-Tf-'  (6) 

where  R  =  AB p^ln{K^I Ki) j p^.  Then,  (4)  becomes 
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D 


dt 


(7) 


^  m  tib-' 


2i 


D  = 


Tt 

RT]{-T)B-i 


(8) 


where  T,  is  dT/dx.  Application  of  (8)  requires  a  knowledge  of  the  variation  of  6^  and 
R  with  temperature. 

If  the  water  in  the  thawed  active  layer  or  the  unfrozen  water  in  the  partially 
frozen  active  layer  or  permafrost  moves  with  a  constant  velocity,  v,  then  a  term  of  the 
form  Cv,vdT/dx  (e.g.,  Guymon,  1984),  where  Cw  it  the  volumetric  heat  capacity  of 
water,  would  be  subtracted  from  the  left  hand  side  of  (4).  This  term,  divided  by  C^, 
would  then  enter  into  the  numerator  on  the  right  hand  side  of  (8)  and  then  both  v 
and  Ca  have  to  be  known  in  order  to  evaluate  D.  Until  present,  only  (2)  has  been 
used  to  evaluate  D;  however,  when  unfrozen  water  is  present,  K  and  Ca  are  strongly 
dependent  on  and  (8)  is  more  suitable. 

The  application  of  (2)  for  evaluating  D  requires  computation  of  Tt  and  Tsx  by 
numerical  methods.  Two  models  will  be  discussed. 


Model  I  The  numerical  approximations  of  the  derivatives  Tt  znd  Txx  axe  (Gerald 
and  Wheatley,  1989) 


and 


7»>+l  T*)— 1 

(9) 

'pi  07’2  j.  Ti 

■■‘(AxV 

(10) 

where  the  integers  t  and  j  reference  positions  and  times  for  the  node  of  interest  and 
Ax  and  At  represent  increments  of  depth  and  time  in  the  observation  mesh.  Equation 
(9)  can  also  be  approximated  by  (Carnahan  et.  al.,  1969) 


7»>+l  _  rjij 

r,  =  ''  +  o((A()l.  (11) 

Equations  (9)  and  (11)  are  alternate  expressions  which  can  be  used  for  numerical 
calculations  where  (9)  involves  two  time  steps  and  (11)  involves  only  one  time  step. 
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Let  A,T  =  7y+‘  -  iy-‘,  A„T  =  r/_i  -  2iy  +  r;Vi  “d  A,r  =  r/*^*  - then, 
neglecting  the  terms  0[(Ax)^]  and  0[(At)^],  substituting  (9)  and  (10)  into  (2)  and 
rearranging  yields 


Similarly,  substituting  (10)  and  (11)  into  (2)  yields 


(12) 


D  = 


(13) 


The  application  of  (12)  requires  three  temperature  profiles  with  two  time  steps 
and  (13)  requires  two  temperature  profiles  with  only  one  time  step.  Equation  (12)  has 
been  used  by  McGaw,  et  al.  (1978);  Nelson,  et  al.  (1985);  Zhang  (1989);  Hinkel,  et  al. 
(1990)  and  (13)  by  Zhang  (1989)  to  compute  the  apparent  thermal  diffusivity  in  the 
a>>'.tive  layer  and  permafrost. 

Model  II  Numerical  approximations  (tnmcation  errors)  of  derivatives  can  be  im¬ 
proved  by  retaining  higher  order  terms  in  the  interpolating  polynomial  (Gerald  and 
Wheatley,  1989).  A  Taylor  series  expansion  with  higher  order  terms  was  used  to  obtain 


T,  =  -  8r/-‘  +  «T!*'  -  t;*^)  +  0|(A()*)  (u) 


and 


T..  =  + 161?-,  -  3or/  +  -  r/.,)  +  o|(Ai)‘)  (is) 

Let  AJT  =  T/"*  -  8r/"‘  -I-  8r/'*'‘  -  and  A'„r  =  -T/.j  +  16r/_i  -  30r,-'  -l- 
1  ~~  neglecting  higher  order  terms  and  substituting  (14)  and  (15)  into 

(2), 


(Ax)^\  A;r 
At  /A'„r' 


(16) 


While  application  of  (16)  requires  five  temperature  profiles  with  four  time  steps, 
it  should  yield  a  more  accurate  estimate  of  D. 
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In  addition  to  the  obvious  requirement  that  the  heat  flow  must  be  conductive, 
there  are  several  general  considerations  for  applying  (12),  (13),  or  (16)  to  field  data. 
First,  the  heat  conduction  equation  has  been  derived  from  continuum  theory  (Carslaw 
and  Jaeger,  1959)  and  is  valid  for  an  infinitesimally  small  volume  or  layer  with  con¬ 
ductive  heat  flow  while  (12),  (13),  and  (16)  are  for  a  finite  volume.  More  accurate 
results  using  these  numerical  equations  may  be  expected  when  Ax  and  At  are  small; 
however,  both  Ax  and  At  are  restricted  by  requirements  associated  with  the  measur¬ 
ing  apparatus.  Second,  the  soil  volume  in  the  active  layer  must  be  in  either  a  frozen 
or  a  thawed  state  for  the  time  steps  of  the  calculations  and  the  soil  volume  should  be 
homogeneous.  Third,  interpolation  errors  in  the  numerical  procedure  depend  upon  the 
degree  of  the  polynomial  used  to  obtain  the  derivatives  so  that  (16)  may  be  expected 
to  give  more  accurate  results  than  (12)  which  is  more  accurate  than  (13).  Round  off 
errors  increase  as  Ax  and  At  are  decreased  since  this  requires  adding  and  subtracting 
temperatures  that  are  more  nearly  the  same  value.  Generally,  round  off  errors  can  be 
reduced  by  using  double  precision  (Gerald  and  Wheatley,  1989).  Fourth,  the  accuracy 
of  the  field  temperature  measurements  influences  the  estimate  of  D.  The  estimate  of 
A(T,  AfT  and  A'(T  are  nearly  independent  of  temperature  accuracy,  since  only  the 
temperature  change  at  one  sensor  is  required,  while  the  estimate  of  AtxT  and  A'^^T 
depends  upon  the  accuracy  of  the  sensors.  The  quantities  AtxT  and  at  any 

depth  should  be  greater  (preferably  much  greater)  than  errors,  6{AxxT)  and  ^(A^^T), 
in  their  determination.  Thus,  temperature  variations  at  depth  should  be  sufficient  to 
make  A^T  >>  ^(A**r)  or  A'^^T  >>  ^(A'„r).  Table  1  shows  6{AxxT)  and  ^(Aj,r) 
for  errors,  6T,  in  the  temperature  measurements. 


Table  1:  Effects  of  errors,  6T,  in  temperature  measurements  on  ^(AzzT)  and 
S{A'„T)  («C) 


6T 

0.2 

0.1 

0.05 

0.01 

S{A„T) 

0.5 

0.3 

0.13 

0.05 

0.03 

7.5 

3.8 

1.88 

0.75 

0.38 

Since  model  II  uses  higher  order  terms  to  calculate  the  derivatives  in  the  expression 
for  D,  it  would  appear  to  be  superior  to  model  I.  An  error  analysis  based  upon  the 
accuracy  of  the  temperature  measurements  required  to  evaluate  (12)  and  (16)  (Table  2) 
shows  that  model  I  produces  somewhat  more  accurate  values  for  D.  For  the  conditions 
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in  Table  2,  the  accuracy  of  the  temperature  measurements  must  approach  ±0.01°C 
or  better  to  obtain  reasonable  estimates  of  D.  Considering  the  increased  complexity 
and  reduced  accuracy  of  model  II,  because  of  measurement  errors,  it  may  not  seem 
worthwhile  to  use  it.  However,  the  reduction  in  truncation  errors  for  T(  emd  in 
model  II  is  a  much  greater  factor  than  the  measurement  errors  which  ultimately  makes 
model  II  much  more  precise  than  model  I.  It  will  also  be  shown  that  the  use  of  model  II 
helps  to  clarify  the  interpretation  of  results  obtained  with  model  I.  In  addition,  model 
II  can  be  used  in  some  cases  when  model  I  becomes  unstable  and  produces  spikes  in 
the  values  for  D. 


Table  2:  Measurement  errors  in  calculated  values  of  D  using  models  I  and  II  at 
a  depth  of  0.3  m  for  a  daily  temperature  wave  with  a  surface  amplitude  of  4°C  with 
Ax  =  0.1  m  and  At  =  10  min. 


Temperature 

accuracy 

(”C) 

Model  I 

Error  in  D 
(mVyr) 

Model  II 

Error  in  D 
(mVyr) 

0.1 

78 

105 

0.05 

39 

52.5 

0.01 

7.9 

10.5 

0.001 

0.9 

1.2 

Synthetic  temperature  time  series  were  calculated  for  surface  temperatures  con¬ 
sisting  of  daily  waves,  annual  waves  with  a  superimposed  daily  variation  zmd  annual 
waves  with  a  superimposed  multi-year  linear  trend.  Wave  parameters  were  chosen  to 
prevent  the  possibility  of  phase  change.  An  input  value  of  D=25  yr~’  was  used 
in  all  of  these  calculations.  Calculated  values  of  D  from  these  synthetic  temperature 
profiles,  using  model  I  (Eq.  12)  and  model  II,  were  then  compared  to  the  input  value 
of  D.  The  stated  errors  are  the  differences  between  the  calculated  and  input  values 
divided  by  the  input  values  (percent  difference).  The  results  help  to  illustrate  the  gen¬ 
eral  considerations  noted  above  and  suggest  additional  factors  that  should  be  taken 
into  account  when  calculating  D  from  a  measured  temperature  time  series. 

RESULTS  AND  DISCUSSION 

Daily  Surface  Temperature  Variations 

Under  stable  weather  conditions,  the  ground  surface  is  subjected  to  a  daily  cycle 
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of  warming  and  cooling,  which  can  cause  measurable  temperature  changes  to  depths 
of  about  one  to  two  meters.  However,  the  presence  of  a  thaw  or  a  freeze  front  shields 
the  underlying  ground  from  these  changes.  Daily  temperature  variations  can  be  used 
to  determine  D  in  a  thawed  or  frozen  active  layer  and  in  the  permafrost  underlying  a 
frozen  active  layer.  Synthetic  temperatures  for  the  daily  cycle  of  warming  and  cooling 
were  generated  with  a  ground  surface  temperature 


r(0,  t)  =  Td  +  AdCos{uJdt)  (17) 

where  Td  is  the  meam  daily  surface  temperature.  Ad  is  the  amplitude  of  the  daily 
surface  temperature,  Ud  =  2v/Pd,  and  the  period,  Pd  =  ^  day.  It  was  assumed  that 
Td  =  -10“C  and  Ad  =  4“C.  Temperatures  were  generated  at  0.1  m  increments  from 
the  surface  to  the  0.5  m  depth  at  10  minute  intervals  using  the  steady-state  solution 
of  (2)  (Carslaw  amd  Jaeger,  1959).  Figure  1  shows  temperature  variations  at  depths 
of  0.2  m,  0.3  m  amd  0.4  m,  the  vairiations  in  AtT  and  AxzT,  and  the  calculated  values 
(using  models  I  and  II)  for  D  at  the  0.3  m  depth  with  At  =  10  minutes.  For  model  I, 
the  calculations  were  made  using  the  temperatures  at  0.2  m,  0.3  m  and  0.4  m  depths, 
while  for  model  II  the  temperatures  from  0.1  m  to  0.5  m  were  used. 

Figure  1(G)  shows  that  there  are  ngnificant  errors  in  values  for  D  calculated  with 
model  I.  The  lairgest  errors  occur  near  times  when  the  temperature  vs.  time  curve  at 
the  0.3  m  depth  passes  through  a  maximum  or  minimum  (i.e.,  when  the  magnitude  of 
A(T  and  AxxT  become  small).  Lau^e  positive  and  negative  spikes  in  D  aire  produced 
at  these  times.  The  best  agreement  of  D  with  the  input  value  occurs  (neav)  where  the 
magnitudes  of  A(T  and  A^sT  aire  largest.  Calculations  of  D  at  other  depths  show  that 
the  spikes  in  D  change  with  depth  and  time  due  to  the  time  lag  as  the  temperature 
wave  penetrates  into  the  ground  &om  the  surface. 

The  results  of  calculations  using  model  II  (Fig.  1C)  aire  much  improved  over 
those  of  model  I.  Spikes  have  been  eliminated  and  the  maximum  error  is  only  4%.  This 
reduction  in  error  appears  to  be  due  to  the  more  accurate  higher  order  approximations 
(smaller  truncation  errors)  for  the  derivatives  in  (2).  Results  for  calculations  using 
double  precision  (sixteen  significant  figures)  were  almost  identical  to  calculations  using 
single  precision  suggesting  that  round  off  errors  were  small  for  both  models  I  and  II 
for  the  conditions  of  these  calculations.  Practical  disadvantages  of  model  II  are  that 
more  data  are  required  and  these  data  need  to  be  very  accurate. 

Figure  2  shows  some  of  the  effects  of  truncation  errors  in  D  when  models  I  and 
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II  are  used  with  difTerent  choices  for  Ax  and  At.  In  Fig.  2B  and  2C,  At  is  constant 
(1  min)  and  Ax  changes  &om  0.1  m  to  0.01  m.  For  model  I,  the  errors  in  D  decrease 
at  all  times  for  this  change  in  Ax  including  a  decrease  in  the  magnitude  and  duration 
of  the  spikes.  In  Fig.  2A  and  2C,  Ax  is  constant  (0.01  m)  and  At  changes  from  10 
min  to  1  min.  Errors  in  D  for  this  change  in  At  decrease  as  At  decreases  except  at 
spikes.  The  duration  of  the  spikes  is  the  same  while  their  magnitude  increases  as  At 
decreases.  The  results  for  model  II  are  similar  but  the  errors  in  D  sire  much  smaller 
(see  Fig.  2B)  and  the  spikes  are  almost  eliminated.  Since  the  conditions  were  the  same 
for  both  models,  this  suggests  that  the  above  errors  are  primarily  truncation  errors  in 
AtT  and  A^^T. 


Daily  and  Seasonal  Surface  Temperature  Variations 

The  ground  surface  temperature  also  varies  with  changes  in  the  weather  and 
seasonal  changes  over  the  year.  These  changes  provide  an  opportunity  to  obtain  D  &om 
depths  of  about  1  m  to  20  m  or  more.  However,  daily  surface  temperature  variations 
are  superimposed  on  these  changes  and  must  be  considered  in  the  methodology.  An 
annual  wave  with  a  superimposed  daily  wave  was  used  to  investigate  this  problem 
Synthetic  temperature  profiles  were  generated  using 


r(0, /)  =  Ty  +  A4Cos{u;dt)  +  AyCos{uy  t)  (18) 

where  Ty  is  the  mean  annual  temperature.  The  second  term  on  the  right  hand  side 
of  (18)  is  the  daily  surface  temperature  wave  and  the  third  is  the  annual  surface  wave 
with  amplitude.  Ay,  and  Uy  =  27r/P,  where  Py  =  1  year. 

Several  sets  of  calculations  of  D  were  made  using  (18)  with  Ad  ranging  from  0°C 
to  4®C  and  A/  from  10  min  to  1  day  with  Ay  =  16®C,  Ty  =  -20®C,  and  Ax  =  0.1  m. 
For  both  models  I  and  II  at  depths  up  to  1.2  m  where  daily  temperature  variations  are 
significant,  using  At  =  1  day,  errors  in  D  were  generally  large  (10  •  30%)  and  depended 
on  Ad.  Using  At  =  10  min,  for  depths  up  to  0.5  m,  the  errors  in  D  for  both  models 
were  similar,  as  expected,  to  those  shown  in  Fig.  IC.  These  results  illustrate  the 
requirement  at  these  shallow  depths  that  At  must  be  much  smaller  than  the  period, 
Pd,  and  show  that  the  efiect  of  the  annual  temperature  wave  on  D  is  not  significant. 

Using  model  I  at  a  depth  of  0.8  m  with  At  =  1  hour,  the  errors  in  D  were  small 
except  near  the  maximum  and  minimum  temperatures  at  this  depth  where  multiple 
spikes  were  produced.  These  spikes  are  a  result  of  daily  temperature  variations  which 


9 


cause  both  A^jT  and  AtT  to  oscillate  with  a  frequency  of  one  day.  The  amplitude  of 
the  oscillations  is  such  that  both  A^T  and  A(T  repeatedly  pass  through  zero  values 
which  produces  positive  and  negative  spikes  in  values  for  D.  Both  the  magnitude  and 
duration  of  the  spikes  were  reduced  when  using  model  II. 

Vadues  for  D  were  calculated  using  models  I  and  II  for  depths  from  4  m  to  20 
m  with  =  4°C,  Ay  =  16®C,  Ty  =  — 20®C,  Ax  =  1  m  and  At  =  5  days  (Fig.  3). 
For  model  I,  the  majumum  percent  difference  between  the  calculated  values  and  the 
input  value  in  Fig.  3  was  less  than  20%.  Diffusivity  values  calculated  by  model  II 
were  almost  identical  to  the  input  value.  Additional  calculations  show  that  D  can  be 
obtained  in  the  depth  range  from  1  m  to  20  m  or  more  with  Ax  ranging  from  a  few 
tenths  of  a  meter  to  a  few  meters  and  At  from  one  day  to  a  few  weeks,  depending 
upon  the  depth  of  interest. 

Climatic  Fluctuations  and  Ground  Surface  Temperature  Variations 

Fluctuations  in  climate  on  a  scale  of  years  produce  changes  in  air  temperatures 
and  in  snow  cover  which  can  affect  permafrost  temperatures.  These  effects  can  be 
monitored  at  depths  below  the  seasonal  changes  by  meztsuring  the  temperatures  in  a 
drill  hole  over  several  years.  While  the  temperature  changes  at  these  depths  are  small, 
they  can  usually  be  measured  with  precision  equipment,  and  D  can  be  estimated  using 
temperature  recording  frequencies  from  a  few  months  to  several  years.  The  method  is 
illustrated  using  the  surface  temperature 


T{0,t)  =  Ty  +  AyCOS{ujyt)  +  Cot  (19) 

where  the  initial  mean  ground  surface  temperature  Ty  =  -20°C.  The  second  term  on 
the  right  side  of  (19)  is  the  annual  temperature  wave  and  the  third  term  is  an  assumed 
linear  temperature  trend  at  the  surface,  where  Cq  is  the  rate  of  temperature  change 
at  X  =  0.  The  initial  condition  is  T(x,  0)  =  +  Ggi,  where  the  temperature  gradient 

Go  =  0.019  “C/m. 

Temperature  profiles  after  10  years  were  generated  with  Ay  =  16“C,  Co  = 
0.2°Cyr~',  for  1  m  space  amd  1  year  time  intervals.  Figure  4  shows  the  calculated 
temperature  profiles  and  values  for  D  predicted  by  models  I  and  II  with  Ax  =  4  m 
and  At  =  1  year.  The  errors  in  D  below  35  m  were  less  than  2%.  The  relatively 
large  errors  above  35  m  were  caused  by  the  annual  temperature  variations  and  can  be 
eliminated  by  choosing  At  <  <  1  year.  At  about  64  m,  round  off  errors  begin  to  affect 
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the  results  for  both  models.  Additional  calculations  show  that  relatively  good  results 
can  be  obtained  for  At  ranging  &om  less  than  a  year  to  several  years. 

APPLICATIONS 

The  above  methodology  suggests  that  it  should  be  possible  to  obtain  values  for 
the  thermal  diffusivity  from  a  temperature  time  series,  not  only  in  the  active  layer 
and  near-surface  permafrost,  but  to  any  depth  where  there  are  measurable  temper¬ 
ature  changes.  However,  the  above  results  are  largely  based  on  assumed  periodic 
surface  temperature  variations  with  fixed  amplitudes  in  soil  where  there  is  no  freez¬ 
ing  or  thawing.  Additional  opportunities  for  determining  D  arc  created  by  weather 
patterns  and  by  the  presence  of  a  freezing  or  thawing  front  in  the  active  layer.  For 
example.  Climate  in  Alaska  and  other  permafrost  areas  is  often  donoinated  by  aperi¬ 
odic  temperature  fluctuations  with  varying  amplitudes  lasting  from  a  few  days  to  a 
few  weeks.  Solar  radiation  at  these  high  latitudes  exhibits  strong  seasonal  variations. 
A  seasonal  snow  cover  attenuates  air  temperature  variations  reducing  the  variations 
in  ground  surface  temperatures.  A  freezing  or  thawing  front  in  the  active  layer  shields 
the  underlying  active  layer  and  permafrost  from  ground  surface  temperature  variations 
since  the  temperature  at  a  front  is  generally  constrained  by  the  requirement  for  phase 
equilibrium  to  be  near  0*’C.  As  the  ground  surface  warms  and  the  active  layer  begins 
to  thaw,  the  remuning  frozen  active  layer  and  permafrost  warm  slowly  with  time. 
As  the  ground  surface  cools  and  begins  to  freeze,  the  remaining  thawed  active  layer 
becomes  isothermal  since  it  is  bounded  on  the  top  by  the  freezing  front  and  on  the 
bottom  by  the  permafrost  table.  Upward  freezing  at  the  base  of  the  active  layer  causes 
the  adjacent  permafrost  to  cool  while  the  deeper  permafrost  continues  to  warm.  The 
temperature  changes  brought  about  by  these  thawing  and  freezing  effects  generally 
continue  for  about  six  to  nine  months  each  year. 

Applications  of  the  method  require  the  acquisition  of  very  precise  temperature 
data  (preferably  d:0.01'’C  or  better)  at  appropriate  depth  (Az)  and  time  {At)  in¬ 
tervals.  Selection  of  values  for  Ax  and  At  must  take  into  account  the  accuracy  of 
the  temperature  measurements,  duration  and  amplitude  of  the  temperature  changes, 
depth  where  D  is  to  be  determined,  and  the  diffusivity  of  the  soil.  Since  the  diffusivity 
of  the  soil  can  easily  vary  from  10  to  50  yr~\  a  factor  of  five  or  more,  it  is  desirable 
to  have  information  which  allows  a  rough  initial  estimate  of  D  to  be  made.  A  periodic 
surface  temperature  wave  can  be  helpful  in  making  initial  estimates  for  Ax  and  At 
provided  some  information  is  available  on  the  wave  parameters.  For  a  periodic  sur¬ 
face  temperature  variation  with  surface  amplitude.  A,  and  frequency,  u;,  the  depth,  X, 
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where  the  amplitude  of  the  temperature  variation  is  Ax  can  be  obtained  from  (Carslaw 
and  Jaeger,  1959) 


Table  3  gives  the  depths  for  a  daily  temperature  wave  {Xd)  and  an  annual  wave 
(Xg)  where  the  amplitude  has  been  reduced  to  Ax  for  D  =  25  Ad  =  4°C, 

and  Ag  =  16°C. 


Table  3:  Depths  for  a  daily  temperature  wave  (Xd)  and  an  annual  wave  (Xy) 
where  the  amplitude  has  been  reduced  to  Ax. 


A.  CO 

0.200 

0.100 

0.050 

0.020 

0.010 

pool 

Xd{m) 

0.34 

0.44 

0.54 

0.68 

0.78 

..12 

Xg  (m) 

12.36 

14.32 

16.27 

18.86 

20.81 

27.31 

Experience  with  models  I  and  II  suggests  that,  in  general.  Ax  and  At  should  be 
made  small  with  Ax  <<  X,  the  depth  of  interest.  At  <<  P  (period  or  duration  of 
sxurface  temperature  change),  A,*r  >>  ^(A„T)  and  A',,T  >>  5(A'„r).  For  the 
active  layer,  these  considerations  and  Table  3  suggest  that  Ax  can  range  from  a  few 
centimeters  to  a  tenth  of  a  meter  or  so  and  At  can  range  from  a  few  minutes  to  an 
hour  or  so.  At  deeper  depths  but  within  the  depth  of  annual  temperature  variations. 
Ax  can  range  from  a  few  tenths  of  a  meter  to  about  one  meter  and  At  from  one  day 
to  several  weeks. 

SUMMARY 

The  purpose  of  this  paper  is  to  investigate  numerical  methods  for  determining  the 
in  situ  thermal  diffusivity,  D,  of  the  active  layer  and  permafrost  from  a  temperature 
time  series.  When  the  active  layer  and  permafrost  contain  unfrozen  water,  the  thermal 
parameters  become  strongly  dependent  upon  temperature.  In  this  case,  it  is  shown 
that 


£>  = 


T, 

Txx  -  PI^(-T)»-> 


(8) 


where  T  is  the  temperature,  Tt,  Tx  and  Txx  the  numerical  approximations  of 
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dTjdt^  dTfdxy  and  d’^Tjdx'^^  respectively,  B  is  a  constant  and  R  is  defined  in  the 
text.  Application  of  (8)  requires  a  knowledge  of  the  variation  of  unfrozen  water, 
and  R  with  temperature. 

Previous  investigators  have  used  the  first  order  numerical  expression  for  D  for  the 
case  of  constant  thermal  parameters,  termed  model  I,  which  is 


where  Ax  is  the  space  step,  At  is  the  time  step,  AjT  =  T/'*'*  -  A„r  = 

r/_,  -  2Tl  +  and  the  integers  i  and  j  reference  positions  and  time.  A  second 
order  approximation  of  the  derivatives  T(  and  Tgg  was  developed  using  a  Tkylor  series 
expansion,  termed  model  II,  which  yields 


(Ax)^\  A;r 
At  yA'„r 


(16) 


where  A;T  =  T/ 
~  ^/+2* 


r-2 


-8T/ 


r-» 


+  8r/+*  -  and  A'„r  =  -  307/  + 


In  addition  to  the  obvious  requirement  for  conductive  hes.t  flow,  they  are  sev¬ 
eral  considerations  in  applying  (12)  and  (16)  to  determine  D.  The  spatial  separation. 
Ax,  of  the  temperature  measurements  (sensors)  and  the  time  intervals.  At,  between 
measurements  must  be  relatively  small  and  the  soil  in  either  a  lirozen  or  thawed  state. 
Estimates  of  AtT  and  A', 7  are  nearly  independent  of  sensor  accuracy  while  Ass7  and 
A',, 7  depend  on  sensor  accuracy.  Temperature  changes  at  depth  should  be  sufficient 
to  make  A(s7  and  A',,7  much  greater  than  the  errors  in  their  determination. 


The  higher  order  approximations  used  for  7|  and  Tgg  in  model  II  result  in  smaller 
truncation  errors  which  makes  model  II  much  more  accurate  than  model  I.  However, 
there  are  also  measurement  errors  associated  with  the  temperatures,  positions,  and 
time.  An  error  analyns  shows  that  the  measurement  errors  are  dominated  by  errors 
in  the  temperature  measurements.  For  a  given  level  of  accuracy  in  the  temperature 
measurements,  the  measurement  errors  for  model  I  are  somewhat  less  than  those  for 
model  II.  The  analysis  also  shows  that  the  accuracy  of  the  temperature  measurements 
must  approach  ±0.01  "C  to  obtain  reasonable  estimates  of  D. 


Additional  information  on  the  application  of  model  I  and  model  II  was  obtained 
by  using  surface  temperatures  consisting  of  a  daily  wave,  an  annual  wave  with  a  super- 
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imposed  daily  wave,  and  an  annual  wave  with  a  superimposed  lineair  trend  to  generate 
synthetic  temperature  time  series  in  the  ground.  Calculated  values  of  D  using  (12) 
and  (16)  were  then  compared  to  the  input  value  used  to  calculate  the  synthetic  tem¬ 
perature  profiles.  The  percent  difference  between  the  calculated  values  and  the  input 
value  is  the  error  in  the  calculations. 

For  a  daily  wave,  model  I  produces  large  spikes  in  the  calculated  values  for  D. 
These  spikes  occur  near  times  where  T  is  a  maximum  or  a  minimum  and  A(T  and 
Azz^  become  very  small.  Using  model  II,  the  spikes  disappear  and  the  errors  are  also 
much  smaller  at  other  times.  These  improved  estimates  using  model  II  are  due  to  the 
reduction  in  truncation  errors  compared  to  model  I.  Round  off  errors  do  not  appear 
to  be  significant  for  either  model  for  the  parameters  used  in  the  calculations. 

For  a  daily  wave  superimposed  on  an  annual  wave,  the  effect  of  the  annual  wave 
is  not  significant  at  shallow  depths  (up  to  0.5  m  in  the  example  calculation).  At  these 
depths,  the  results  for  D  were  similar  to  those  for  a  daily  wave  only.  As  the  amplitude 
of  the  daily  wave  increased,  the  errors  in  D  increased  and  as  the  depth  increased  to  1.2 
m  and  beyond  the  errors  decreased.  At  a  depth  of  0.8  m,  calculated  values  for  D  using 
model  I  resulted  in  a  series  of  spikes  near  the  temperature  maximum  and  minimum 
and  normal  values  elsewhere.  These  spikes  were  produced  by  the  daily  temperature 
variations  which  caused  the  temperature  dependent  quantities  on  the  right  sides  of 
(12)  and  (16)  to  oscillate  about  zero  for  a  period  of  time.  The  number  of  spikes,  their 
amplitude,  and  duration  were  substantially  reduced  using  model  II.  At  a  depth  of 
several  meters  to  20  m  or  more,  model  I  produced  values  for  D  which  were  in  error 
by  less  than  20%.  However,  when  model  II  was  used,  the  values  for  D  were  nearly 
identical  to  the  input  value. 

An  annual  surface  wave  with  a  superimposed  linear  trend,  which  could  be  caused 
by  climatic  fluctuations,  was  used  to  show  that  D  could  be  obtained  at  depths  below 
those  of  annual  variations  (about  35  m  for  this  example).  The  errors  in  values  for  both 
modeb  I  and  II  were  less  than  2%  with  model  II  being  more  precise. 

Values  for  D  can  be  obtained  at  any  depth  where  measurable  temperature  changes 
occur.  Short  term  weather  patterns,  the  presence  of  a  freezing  or  thawing  front  in  the 
active  layer,  and  longer  term  climatic  fluctuations  and  trends  create  opportunities 
for  determining  D.  Application  of  the  method  requir'e*  the  acquisition  of  very  precise 
temperature  data  (preferably  ±0.01*’C  or  better)  at  appropriate  depth  (Ax)  and  time 
(At)  intervals.  Selection  of  values  for  Ax  and  At  must  take  into  account  the  accuracy 
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of  the  temperature  measurements,  duration  zmd  amplitude  of  the  temperature  changes, 
depth  where  D  is  to  be  determined,  and  the  diifusivity  of  the  soil.  In  general.  Ax  <<  X 
(the  depth  of  interest),  At  <<  P  (  period  or  duration  of  surface  temperature  changes), 
and  AixT  and  A'^^T  much  greater  than  the  errors  in  their  determination.  For  the 
active  layer,  Ax  can  range  from  a  few  centimeters  to  a  tenth  of  a  meter  or  so  and 
At  from  a  few  minutes  to  an  hour  or  so.  At  deeper  depths  in  the  perm2ifrost  but 
within  the  depth  of  annual  temperature  variations.  Ax  can  range  &om  a  few  tenths  of 
a  meter  to  a  meter  or  more  and  At  from  one  day  to  several  weeks.  Below  the  depth 
of  annual  temperature  variations,  values  for  Ax  of  several  meters  and  At  of  a  year 
or  more  may  be  appropriate.  The  method  is  not  limited  to  permafrost  regions  and 
should  be  applicable  to  any  soil  or  rock  where  the  heat  flow  is  conductive. 
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Figure  Captions 


Figure  1.  (A)  Synthetic  temperature  time  series  at  three  depths  for  a  daily  siirface 
temperature  wave  with  a  surface  amplitude  of  4'’C,  (B)  and  AtT  at  the 

central  depth  of  0.3  m,  and  (C)  a  comparison  of  the  calculated  values  of  the 
thermal  diffusivity  at  the  0.3  m  depth  with  the  input  value  of  25  t/r~*.  For 

these  calculations,  Ax  =  0.1  m,  At  =  10  min. 

Figure  2.  Variations  of  the  calculated  thermal  diffusivity  at  a  depth  of  0.3  m  for  different 
values  of  Ax  and  At:  (A)  Ax  =  0.01  m  and  At  =  60  min,  (B)  Ax  =  0.10  m  and 
At  =  1  min,  and  (C)  Ax  =  0.01  m  and  At  =  1  min. 

Figure  3.  (A)  Five  synthetic  temperature  profiles  for  depths  below  the  daily  zero 
temperature  variation  and,  (B),  a  comparison  of  the  calculated  values  between 
models  I  and  II.  For  these  calcidations.  Ad  =  4°C,  Ax  =  1  m,  and  At  =  5  days. 

Figure  4.  (A)  Five  synthetic  temperature  profiles,  czdculated  after  ten  years,  for  a 
warming  trend  of  ground  surface  temperature  with  a  superimposed  annual  surface 
temperature  wave  and,  (B),  a  comparison  of  the  calculated  thermal  diffusivities 
&om  the  temperature  profiles.  For  the  temperature  calculations.  Ay  =  16°C, 
Ax  =  1  m,  and  At  =  1  year,  and  for  the  thermal  diffusivity  calculations.  Ax  =  4 
m  and  At  =  1 


year. 


0 


400 


LEGEND 
Modrt  I 
Modol  II 


800 

(min) 


1000 


DEPTH  (m) 


TEMPERATURE  (C)  DIFBJSIVITY  (m*2/yr) 


DEPTH  (m) 


